function spheresurf(data,scale)
%SPHERESURF Display a matrix of surface data overlaid on a sphere.
%   SPHERESURF(DATA,SCALE) displays the heatmap, elevation, or antenna 
%   response pattern surface DATA, wrapped onto the surface of a sphere,
%   with height SCALE relative to the sphere's radius.  The input DATA is
%   typically a matrix of values expressed over a Cartesian grid of 
%   latitude and longitude angles.  SCALE must be a positive number or 
%   zero.  For SCALE = 0, the data is shown as a flat heatmap on the 
%   surface of the sphere.
%
%   SPHERESURF(DATA) displays the data with relative magnitude SCALE = 1.
%
%   Example 1: Display the 3-D response pattern generated by a 6-element 
%   uniform rectangular array.
%
%       load('arrayresponse.mat')
%       hf = figure('Visible','off','Position',[400 260 1250 430]);
%       axes('Units','pixels','Position',[50 60 350 350])
%       spheresurf(r,0)
%       title('Scale = 0')
%       axes('Units','pixels','Position',[450 60 350 350])
%       spheresurf(r,0.3)
%       title('Scale = 0.3')
%       axes('Units','pixels','Position',[850 60 350 350])
%       spheresurf(r,inf)
%       title('Scale = inf')
%       movegui(hf,'center')
%       set(hf,'Visible','on')
%
%   Example 2: Display the surface of the Earth with exaggerated elevation
%   values.
% 
%       load('topo.mat','topo','topomap1')
%       figure
%       spheresurf(topo,0.1)
%       colormap(topomap1)

% Paul Fricker 11 May 2012 The MathWorks, Inc.

% Input checking
if nargin == 0
    error('')
elseif nargin == 1
    scale = 1;
elseif nargin > 2
    error('respsurf:nArgs','Incorrect number of input arguments')
end

if ~isnumeric(data)
    error('respsurf:isDataNumeric','The first input (data) must be a numeric array.')
end
if ~isreal(data)
    error('respsurf:isDataReal','The input data must be real.')
end
if ~ismatrix(data)
    error('respsurf:isDataMatrix','The first input (data) must be a matrix.')
end
if any(size(data)==1)
    error('respsurf:isDataVector','The first input (data) must be a matrix.')
end

if ~isnumeric(scale)
    error('respsurf:isScaleNumeric','The second input (scale factor) must be numeric.')
end
if ~isreal(scale)
    error('respsurf:isScaleReal','Scale factor must be a real number.')
end
if scale < 0
    error('respsurf:isScalePositive','Scale factor must be a positive number (or zero).')
end

% Create a unit sphere
[x,y,z] = sphere(size(data));

% Scale the data and merge it with the sphere
if ~isinf(scale)
    r = (1 + scale*data/max(abs(data(:))))/(1 + scale);
else
    r = data/max(abs(data(:)));
end
xmod = r.*x;
ymod = r.*y;
zmod = r.*z;

if any(r(:) <= -1/(1+scale))
    warning('respsurf:negativeWarning', ...
            ['The provided input DATA and SCALE factor produce negative ' ...
             'values which are visible through the back of the sphere.'])
end

% Create the surface plot
ha = newplot;
surf(xmod,ymod,zmod,data,'Parent',ha), shading interp
set(ha,'XTickLabel','','YTickLabel','','ZTickLabel','', ...
       'XLim',[-1 1],'YLim',[-1 1],'ZLim',[-1 1], ...
       'TickLength',[0 0])
axis square
colorbar


function [x,y,z] = sphere(nm)
%SPHERE Create a unit sphere.
%   [X,Y,Z] = SPHERE(NM) creates an N-by-M set of (X,Y,Z)-values on the
%   unit sphere.

n = nm(1);
m = nm(2);

% Build vectors of angles
theta = (-(m-1):2:(m-1))/(m-1)*pi;
phi   = (-(n-1):2:(n-1))'/(n-1)*pi/2;

% Compute (X,Y,Z)-coordinate locations
sintheta = sin(theta); sintheta([1 m]) = 0;
cosphi   = cos(phi);   cosphi([1 n])   = 0;
x = cosphi*cos(theta);
y = cosphi*sintheta;
z = sin(phi)*ones(1,m);
